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Abstract 

For the stationary advection-diffusion problem the standard continuous Galerkin method is un- 
stable without some additional control on the mesh or method. The interior penalty discontinuous 
Galerkin method is stable but at the expense of an increased number of degrees of freedom. The 
hybrid method proposed in [5] combines the computational complexity of the continuous method 
with the stability of the discontinuous method without a significant increase in degrees of freedom. 
We discuss the implementation of this method using the finite element library deal . ii and present 
some numerical experiments. 

1 Introduction 

We consider the advection-diffusion equation 

(1.1) -eAu + b- Vu = f mncR d 

(1.2) u = g on dfl 

with < s <C 1, b € W°°(div, O), / G L 2 (f2) and g £ H^ 2 (p,). For simplicity we assume the region Q 
is polygonal. We also assume p := — iV- b > and then we have a weak solution u G H 1 ^!). 

It is well known that this problem can exhibit boundary or internal layers in the convection 
dominated regime and that for the standard continuous Galerkin (cG) formulation these layers cause 
non-physical oscillations in the numerical solution. Several adaptations to the cG method are effective 
but space does not allow their discussion here. We refer readers to [9] for a full description of these 
approaches. Discontinuous Galerkin (dG) methods also offer a stable approach for approximating this 
problem. However the number of degrees of freedom required for dG methods is in general considerably 
larger than for cG methods. 

We describe an alternative approach also studied in [5l El [7] . A dG method is applied on the layers 
and a cG method away from the layers. We call this approach the continuous-discontinuous Galerkin 
(cdG) method. The hypothesis is that provided the layers are entirely contained in the dG region the 
instability they cause will not propagate to the cG region. Note that in our formulation there are no 
transmission conditions at the join of the two regions. 

Here we present the cdG method and discuss its implementation using the deal . ii finite element 
library. We additionally provide some numerical experiments to highlight the performance of the 
method. 
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2 Finite Element Formulation 



Assume that we can identify a decomposition of Q, := Q, cG U Q dG where it is appropriate to apply the 
cG and dG methods respectively. We do not consider specific procedures to achieve this here, but 
generally it will be that we wish all boundary and internal layers to be within f2 dG . Identifying these 
regions can be done a priori in some cases or a posteriori based on the solution of a dG finite element 
method. Consider a triangulation Th of $7 which is split into two regions T h cG and 7^ dG where we will 
apply the cG and dG methods respectively. For simplicity we assume that the regions T h cG and 7^ dG 
are aligned with the regions Q and Q dG and the set J contains edges which lie in the intersection 
of the two regions. Call the mesh skeleton £h and the internal skeleton £?. Define T as the union of 
boundary edges and the inflow and outflow boundaries by 

T in ={x £dQ:b-n<0} 
T out ={x G dQ : b- n > 0} 

where n is the outward pointing normal. Define T cG (resp. T dG ) to be the intersection of V with 7^ G 
(resp. T d ). By convention we say that the edges of J are part of the discontinuous skeleton £ dG and 
:= £h\ £h G - With this convention there is potentially a discontinuity of the numerical solution 
at J. Elements of the mesh are denoted E, edges (resp. faces in 3d) by e and denote by He and h e 
the diameter of an element and an edge, defined in the usual way. 

The jump [ • J and average § • § of a scalar or vector function on the edges in £h are defined as in, 
e.g., [J. 

Definition 2.1. Define the cdG space to be 



(2.2) 



V, 



cdG 



{v G L\n) : v\ E G F k , v\ dQndQcG = g, v\ QcG G H X #T G )} 



where F k is the space of polynomials of degree at most k supported on E. This is equivalent to applying 
the usual cG space on £l cG and a dG space on Q, dG . 

We may now define the interior penalty cdG method: Find Uh £ KdG sucn that for all Vh G VcdG 

B £ (u h ,v h ) = B d (u h ,v h ) + B a (u h ,v h ) = L £ (f,g;v h ) 



where 



Bd(u h ,v h ) 



B a (u h ,v h ) 



^2 / eVu h -Vv h - / (b-Vv h )u h - / (V-b)u h v h 
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^2 / b-lv h }u h + ^2 (b-n)u h v h 



and 



L e (f,g;v h ) 
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Here a is the penalization parameter and d £ {—1,0, 1}. Note that through the definition of V^ G the 
edge terms are zero on ££ G and the method reduces to the standard cG FEM. If we take Th = T dG , 
i.e., the entire triangulation as discontinuous, we get the interior penalty (IP) family of dG FEMs (see, 

eg.', m)- 

The work of [8] shows that the cG method is the limit of the dG method as a — > oo. A reasonable 
hypothesis is that the solution to the cdG method is the limit of the solutions to the dG method 
as the penalty parameter a — > oo on e G ££ G , i.e., super penalising the edges in £f l G . Call a c Q and 



o~dG the penalty parameters for edges in £^ G and £ dG respectively. Call the numerical solution for 
the cdG problem u cc \G,h G KdG- ^ ne solution to the pure dG problem on the same mesh is denoted 
u dG,h G V^ G where V^ G is the usual piecewise discontinuous polynomial space on Th- Then we have 
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Theorem 2.3. The dG solution converges to the cdG solution of (|l.ip - (|1.2p as g c q —¥ oo, i.e., 

lim (ucdGh - u dG) h) = 

We do not prove this result here but direct readers to [I] for a full discussion. Although this result 
does not imply stability of the cdG method (indeed, for the case where the Q cG region is taken to be 
the whole of fi it shows that the cdG method has the same problems as the cG method), but does 
indicate that investigation of the cdG method as an intermediate stage between cG and dG is justified. 
Hence, it aids into building an understanding of the convergence and stability properties of the cdG 
method, based on what is known for cG and dG. This, in turn, is of interest, as the cdG method offers 
substantial reduction in the degrees of freedom of the method compared to dG. 

3 Numerical Implementation 

The cdG method poses several difficulties in implementation. One approach is to use the super 
penalty result of Theorem 12. 3l to get a good approximation to the cdG solution. However this will give 
a method with the same number of degrees of freedom as dG. We therefore present an approach to 
implement the cdG method with the appropriate finite element structure. We discuss this approach 
with particular reference to the deal.ii finite element library [2 [3]. This is an open source C++ 
library designed to streamline the creation of finite element codes and give straightforward access to 
algorithms and data structures. We also present some numerical experiments. 

3.1 Implementation in deal . ii 

The main difficulty in implementing a cdG method in deal, ii is the understandable lack of a native 
cdG element type. In order to assign degrees of freedom to a mesh in deal.ii the code must be 
initialised with a Tri angulation and then instructed to use a particular finite element basis to place 
the degrees of freedom. Although it is possible to initialise a Triangulation with the dG and cG 
regions set via the material_id flag, no appropriate element exists. In the existing deal . ii framework 
it would be difficult to code an element with the appropriate properties. A far more robust approach 
is to use the existing capabilities of the library and therefore allow access to other features of deal . ii . 
For instance without the correct distribution of degrees of freedom the resulting sparsity pattern of the 
finite element matrix would be suboptimal, i.e., containing more entries than required by the theory 
and therefore reducing the benefit of shrinking the number of degrees of freedom relative to a dG 
method. 

The deal.ii library has the capability to handle problems with multiple equations applied to a 
single mesh such as the case of a elastic solid fluid interaction problem. In our case we wish to apply 
different methods to the same equation on different regions of the mesh, which is conceptually the 
same problem in the deal . ii framework. In addition we will use the hp capability of the library. 

The deal . ii library has the capability to create collections of finite elements, hp : : FECollection. 
Here multiple finite elements are grouped into one data structure. As the syntax suggests the usual use 
is for hp refinement to create a set of finite elements of the same type (e.g., scalar Lagrange elements 
FE_Q or discontinuous elements FE_DGQ) of varying degree. Unfortunately it is not sufficient to create 
a hp: : FECollection of cG and dG elements as the interface between the two regions will still be 
undefined. In order to create an admissible collection of finite elements we use FEJJOTHING . This is a 
finite element type in deal . ii with zero degrees of freedom. Using the FESystem class we create two 
vector- valued finite element types (FE_Q, FEJJDTHING ) and (FE_NOTHING , FE_DGQ) and combine them in 
a hp: : FECollection. We apply the first FESystem on the cG region, and the second on the dG region. 
Now when we create a Triangulation initialised with the location of cG and dG elements the degrees 
of freedom can be correctly distributed according to the finite element defined by hp: : FECollection. 

When assembling the matrix for the finite element method we need only be careful that we are 
using the correct element of hp: : FECollection and the correct part of FESystem. The most difficult 
case is on the boundary J where from a dG element we must evaluate the contribution from the 
neighbouring cG element (note that in the cdG method a jump is permissible on J). 
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If we implement the cdG method in deal.ii in this way we create two solutions: one for the 
FE_Q-FE_NDTHING component and another for the FEJJOTHING -FE_DGQ component. Consider a domain 
= (0, l) 2 in R 2 , b = (1, 1) T and i? = —1. The Dirichlet boundary conditions and the forcing function 
/ are chosen so that the analytical solution is 



(3.1) 



u(x, y) = x + y(l — x) + 



This solution exhibits an exponential boundary layer along x = 1 and y = 1 of width 0(e), e = 10~ 6 . 
We solve the finite element problem on a 1024 element grid and fix Q cG = [0, 0.707) 2 . This is larger 
than is required for stability (see Example 1 below) but shows the behaviour of FEJJOTHING more 
clearly. We show each of the components of FE_SYSTEM and the combined solution. For comparison 
we also show the dG finite element solution for the same problem. 



(a) cG-FE_N0THING 



(b) FE_N0THING -dG 



(c) Combined cdG solution 



(d) dG soluti 



Figure 3.1: Solution components of FE_NOTHING implementation applied to Example 1 with e = 10 6 . 

One advantage of following the deal.ii framework is that the data structures will allow the 
implementation of hp methods. In fact we can envisage the implementation of a hpe method where 
at each refinement there is the possibility to change the mesh size, polynomial degree or the element 
type. We propose no specific scheme here but simply remark that implementing a hpe method is 
relatively straightforward with the FE_N0THING approach. 



3.2 Numerical Examples 

We present two numerical experiments highlighting the performance of the cdG method. Both ex- 
amples present layers when e is small enough. In each case we fix the region where the continuous 
method is to be applied then vary e. This causes the layer to steepen. In the advection dominated 
regime, i.e., e large and no steep layer present, we see the cdG solution approximates the true solution 
well. As we make e smaller the layer forms and extends into the continuous region. As e becomes 
smaller still the layer leaves the continuous region and the performance of the dG and cdG method 
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is indistinguishable. In each experiment we pick the £1 and Q regions so that with the given 
refinement the region T dG consists of exactly one layer of elements and coincides with O dG . 

Example 1 Consider again the problem with true solution (|3.ip presented above. We solve the finite 
element problem on a 1024 element grid and fix Q cG = [0, 0.96875) 2 so exactly one row of elements is 
in Q dG . As we vary e the layer sharpens and moves entirely into the dG region. 



0.06 r 0.6 r 




epsilon epsilon 

(a) \\u cdG . h - u dG ,h\\ L 2(n) (b) [|«cdG,h - UdG,Ji[U°°(n) 

Figure 3.2: Decreasing e with a fixed Q decomposition in Example l.The maximum difference in either 
norm occurs when the layer is sharp but not contained entirely in Q dG . 

As we can see from Figure 13.21 before the layer has formed the two methods perform well. As the 
layer begins to form with decreasing e it is not entirely contained in the discontinuous region and the 
error peaks. As the layer sharpens further it is entirely contained in the discontinuous region and the 
difference between the two solutions becomes negligible. 

Example 2 Now we look at a problem with an internal layer. Let the advection coefficient be given 
by b = (—x,y) T and pick the boundary conditions and right hand side / so that the true solution is 

u(x, y) = (1 - y 2 )erf (r^=) > 

where erf is the error function defined by 

erf(x) = [ X e~ t2 dt. 
V 71 " Jo 

We solve on the region Q = (— 1, l) 2 . The solution has an internal layer along y = of width 
0(y/e) and we fix n cG = {(x,y) : x G [-1,-0.0625) U (0.0625, 1], y G [-1,1]}. 

In Figure [331 we notice same the same behaviour as in Example 1. If the layer exists it must be 
contained within the discontinuous region for the two methods to perform equivalently. In Figure 13.31 
we can see the cdG solution for various e with the oscillations clearly visible when e = 10 -4 . When 
the layer is sharpened, the oscillations disappear. 
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(a) e = 1CT 2 (b) e = 1(T 4 (c) e = 1CT 6 

Figure 3.3: The cdG solutions for Example 2 for various e. When e = 10~ 4 the layer is steep enough 
to cause oscillations but not sharp enough to be contained entirely in f] dG . In this case the oscillations 
are clearly visible, but they are not present when e = 10 -6 as the layer has moved entirely within 

n dG . 



0.07 r 0.3 r 
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epsilon epsilon 
(a) ||«cdG,h — WdQ,h||i2(n) (b) ||McdG,/i - UdG,fc||z,°°(n) 

Figure 3.4: Decreasing e with a fixed decomposition in Example 2. As in Figure [3721 the maximum 
difference in either norm occurs when the layer is sharp but not contained entirely in J7 dG . 
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